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Abstract 

We propose a new method based on variational principle for anal¬ 
ysis of photonic crystal (PC) slabs. Most of the methods used today 
treat PC slab as a three-dimensional (3D) crystal and this makes them 
very time and/or memory consuming. In this method we use Bloch 
theorem to expand the field on infinite plane waves which their am¬ 
plitudes depend on the component perpendicular to the slab surface. 

By approximating these amplitudes with appropriate functions, we 
can find modes of PC slabs almost as fast as we can find modes of 
a two-dimensional (2D) crystal. Besides this advantage , we can also 
calculate radiation modes with this method which is not feasible with 
3D Plane Wave Expansion (PWE) method. 

Keywords: variational principle, photonic crystal slab, plane wave expan¬ 
sion 
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1 Introduction 

Since the prediction of forbidden gaps in the band-structure of some peri¬ 
odic dielectric structures that we now call photonic crystals (PCs) [1], these 
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structures have affected many scientific fields and have found many indus¬ 
trial applications. Today, PCs are used to increase solar cells efficiency [2-5], 
improve lasers and light emmiting diodes specs [6-8], and design new kind 
of waveguides and optical hbers [9-13]. These nano-structures are a good 
platform for optical circuits [14, 15] and seem to hnd application in some 
new research helds such as quantum information and quantum computa¬ 
tion [16-21]. 

Periodicity of PCs dielectric constant can be in one, two, or three dimen¬ 
sions [22-24]. Fabrication of full 3D PCs or perfect 2D ones is difficult or 
sometimes impossible. This fact has made scientist to do their best to work 
with slab of PCs because they are fabricated easily with conventional fabri¬ 
cation technologies. PC slabs have periodic permittivity in two dimensions, 
but hnite thickness along the third one. They are sometimes known as hnite 
thickness 2D crystals. These crystals can conhne and guide electromagnetic 
waves according to distributed Bragg reflection in the plane of the slab and 
total internal reflection in the slab normal direction [11,25-28]. Up to now 
many methods have been proposed for analysis of PCs which fall into two 
categories; frequency domain and time domain. Frequency domain methods 
such as hnite element, plane wave expansion, and hnite diherence frequency 
domain deal with phasors of electromagnetic helds, but time domain ones 
like hnite diherence time domain (FDTD) hnd helds evolution according to 
Maxwell’s equation. Since PC slab is not a perfect 2D crystal, we need to 
simulate it as a full 3D one and this has made its analysis time and/or mem¬ 
ory consuming. Some innovative methods have been proposed to solve this 
problem [29-33], but almost all of them have some shortcomings among them 
we can name, lower accuracy and limited frequency range analysis. 

In this paper we hrst introduce our new method which can hnd PC slab 
modes by calculating eigenvalues of a matrix twice as large as the matrix 
appears in 2D PWE method. Then we obtain band structure of a sample 
PC slab and compare it with the results of conventional 3D PWE and FDTD 
methods. We also compare eigen held of this crystal at a high symmetry point 
obtained from the proposed method with that of 3D PWE method. To show 
its ability to analyze PC waveguides we hnd guided modes of a sample PC 
waveguide, and hnally we compare the efficiency of this method against 3D 
PWE method. 
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2 Variational method for PC slab analysis 

According to variational principle one can finds an estimation to smallest 
eigenvalues of a Hermitian operator L, by choosing a suitable trial function 
as the eigen function of the operator and then trying to minimize 

(/|L|/) 

(/I/) ’ 

where |/) is the trial function. In fact the popular PWE method is based on 
the variational principle. In this method the trial function, that is one of the 
electric or magnetic helds, is chosen according to Bloch theorem to be 

F(r) = ^ 

G 

where is a periodic function with the same periodicity as that of the 

dielectric constant of the crystal, k is the Bloch wave vector and G equals 
with bjS being the primitive vectors of the reciprocal lattice. 
The second equality is written by substituting *hK,(i’) with its Fourier series 
expansion. In photonic crystal slab there is no periodicity in the vertical 
direction of the slab surface. Hence to use the standard PWE method we 
need to create an artihcial periodicity in this direction. This means that we 
are estimating the held dependency on the normal component of the slab by 
a Fourier series. Actually this is not a good estimation, because the number 
of Fourier series coefficients that should be determined can be very large. 

Here we want to show we can choose a simpler trial function for the held 
dependency on the normal component of the slab which results in faster 
calculation of crystal modes. For clarity of the formulation written in the 
remaining of the paper, suppose we want to analyze the photonic crystal 
slab shown in hgure 1. This crystal is composed of a triangular lattice of 
circular air columns etched through a dielectric slab with dielectric constant 
er = 11.9, and thickness t = 0.6a, where a is the lattice constant. The 
radius of air columns equals 0.3a. According to Bloch theorem we can write 
magnetic held as 

= ( 1 ) 

where z) is a periodic function of x and y for every value of If we 
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Figure 1: PC slab analyzed in this paper constructed by etching circular air 
columns through a dielectric slab with dielectric constant e,, = 11.9 

substitute ^{r^y^z) by its Fourier series, we can write 

H(r) = 

G 

= (^G||(2;)eG|| + hGx(^)eGx + 

G 

X (2) 

where eGy and Gg^ are unit vectors along and perpendicular to k + G re¬ 
spectively. Here G is the vector in reciprocal lattice of a 2D crystal with the 
same pattern as that of our PC slab. To find a good estimation for hG||(2:), 
hQ,^{z), and hG^{z) we need to survey field distribution in a dielectric slab 
waveguide problem. We know in that problem, amplitudes of guided electric 
and magnetic fields inside the slab are sinusoidal functions of z component. 
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but decay exponentially outside the slab as |; 2 | increases, that is 


F,(r,t) = , (i = ||, ±,;2) 

{ AoSin(A; 2 :) or yleCOs(A; 2 :) , \z\ < t/2 

C^e-»iz-tl 2 ) ^ z>t/2 

Q^^a{z+t/2) ^ ^ ^ 


( 3 ) 


where F stands for electric E, or magnetic H, fields, u is the angular fre¬ 
quency, and Aoi , C„, C;, a, and k are constants. It can be shown 

-I- k-2 

(jj = - and a = K — u , (4) 

where k = |k|. Considering each of the even/odd or TE/TM modes of the 
waveguide we can also write four relations between a and /c. 


j a = /ctan(/ct/2) 

1 ^ 0 ; = —k cot(/ct/2) 

j a = k tan{kt/2)/er 

[ a = —kcot{kt/2)/er 


(5) 


We can simply show that fields prohle do not change outside the slab if we 
carve out the triangular lattice of circular columns, that is electromagnetic 
held decays exponentially outside the PC slab as \z\ increases. If we can 
approximate its prohle inside the PC slab by sinusoidal functions, then our 


5 



estimated trial function in TE-like mode becomes 






{ sm{kG^^z) 

sin(fcG||t/2)e"“G(^“‘/^) 
— sin(A;G|| 


^G||/G||(2:) 

{ sin(fcGx2:) 

sin(A;Gx^/2)e““°*'^“*/^^ 
— sin(/cGxV2)e"®*^^+*/^^ 


kl < t/2 

z > t/2 
z < —t/2 


\z\ < t/2 
2 ; > t/2 
z < —t/2 


— ^G_L/Gi(^) 

hcA^) = -/!«; + GIAgii 

COs(fcG||^)/fcGj| 

sin(fcG t/2) f/n\ 

X / _)_ II g-aG(^-v2) 

' «G 

Qq 

= ^G||/g,(^)- 


kl < t/2 

^ > t/2 
2 ; < —t/2 


( 6 ) 


/tG||(^) and /tGx(^) are chosen to be continuous at \z\ = t/2 and is 

written such that V • H(r) = 0. TM-like mode trial functions can be written 
in the same manner. 

^G||; ^G^; ^G|p ^Gj_; and gg are parameters which can be determined 
by minimizing 

(HWILhIHW) 

{H(r)|H(r)) ’ ^ ' 

where 

^"‘'0 ■ 

This is again a time consuming problem. To further simplify it, suppose we 
have suitable values for og, ^G\\i and Then, minimization of expression 
(7) becomes minimizing (Ag| M |Ag) provided that (Ag| 1^ I^g) = 1, where 


I^g) = 

and N = 



( nil 

oy 

VO 

^ 22 /_ 


[ niii 

mi2\ 

Vg'-g 

\m2i 

^22/_ 


G,G' 


( 8 ) 
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In (8) r^G are Fourier series coefficients of T]{r) = l/er(r), and elements of 
matrices M and N equals 

mu = ecii -ec '(|k + GP(/g/(2)|/g.(2;)) 

+ \^ + G\^\^ + Gf{fodz)\fGM) 

+ (/^(^)I/g]|(^)) + I'^ + G'|^(/g^(^)|/g||(2^))), 

mi2 = bg^ ■ eGji ((/gJi (^)I/gx(^)) 

+ |k + G'p(/G^(2;)|/G||(2:))j, 

m2i = eG|| ■eG'^((/GY(^)l/G||(^)) 

+ |k + {f(.,^{z)\fG^{z))^ , 

17122 = \k + G\\k + G'\{fG'^{z)\fGA^)) 

+ GGx ■ eGY(/GY(^)l/Gi(^))> 

( (/G||(^)|/G||(^)) 

nu={ +A + G\AfGAA\fG.{z)) , G = G' , 

[ 0 , G ^ G' 

„ _ r ifGAAlfGAA) , G = G' 
rr22 - I 0 , G ^ G' • 

Using Lagrange multipliers method we can determine |^g) elements by 
solving the generalized eigenvalue problem, 


M|LIg) + AN|Ag) = 0. (10) 

As said above we hrst need to set the values of og, ^Gy) and A;gj_ 
order to obtain the simplihed eigenvalue problem of (10). For this purpose 
we temporarily assume angular frequency u, to be that of dielectric slab 
waveguide with the same slab thickness and permittivity equals the effective 
permittivity of PC slab. We mean by effective permittivity, the coefficient of 
Fourier series of er(r) with G = 0. Please notice that we only consider the 
hrst band of dispersion diagram of slab waveguide, that is we set 


UJ = wid^l). 

After that using (4) we can write 

ttG = a/|k + Gp - a;2. 
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and then ^Gh and are obtained by solving the following eqnations 


TE-like: 

TM-like: 


ttG = ^Gii tan(fcG||^/2) 
og = -ka^ coi{kGj,/2)/er^^ 
ttG = -^G|| cot(A;G||f/2) 
ttG = ^Gx tan(fcGxf/2)/er^g 


where is the effective permittivity of PC slab. 


3 Results comparison 

The calcnlated TM-like and TE-like band strnctnres of the crystal shown in 
hgnre 1 by the proposed variational method have been compared with the 
resnlts of standard 3D PWE method in hgnres 2 and 3 respectively. For 
PWE method, nnmber of Fonrier series terms in xy plane was bonnded to 
Nx = Ny = 5 and to = 6 in z direction. We also used Nx = Ny = 5 

for variational method. As is seen results have good agreements. In fact 
in TM-like mode our method gives more precise result. Discrepancy in high 
frequency bands of TE-like mode is mostly due to inaccuracy of 3D PWE 
method near the light cone. As mentioned, in 3D PWE method we need to 
make an artihcial periodicity along the slab normal direction by introducing 
inhnite slabs similar to the main one above and below it. This is done because 
held decay exponentially outside the slab and we can assume other slabs do 
not change the held distribution near and inside the main slab. But when 
we go near to the light cone that is when we are on the edge of conhned and 
radiation modes of the slab, the held decay slowly outside the slab and this 
assumption is not valid any more so 3D PWE method loses its accuracy. 

One of the advantages of the proposed method over the 3D PWE one 
is that besides calculating guided modes of PC slab, we can also obtain 
radiation modes by this method, hgnre 4 shows guided and radiation TM- 
like modes of the crystal of hgnre 1 calculated by the variational method and 
compares it with FDTD method result. 

We have also used our method to plot TM-like magnetic held distribution 
in the mid-plane of the PC slab and in a vertical plane to it at the high- 
symmetry point of the reciprocal lattice M*^^\ hfth mode at the M point, and 
compared it with the result of PWE method in hgnre 5. Again we see good 
agreement between the results. 


frequency (a/A) 



Figure 2: Comparison of TM-like band structures of the crystal shown in 
hgure 1 obtained from the variational method by that of 3D PWE method. 
Green region shows the light cone. 


To examine our method capability in dealing with large unit cell crystals, 
we plotted in hgure 6 dispersion diagram of a PC slab waveguide by both 
FDTD and variational methods. This waveguide is built by hlling one row 
of air columns with dielectric in the crystal of hgure 1 in F — K direction. 
It is shown in the upper inset of hgure 6 and has three guided modes which 
our method has found all of them. The lower inset of hgure 6 shows the 
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Figure 3: Comparison of TE-like band structures of the crystal shown in 
hgure 1 obtained from the variational method by that of 3D PWE method 


unit cell used to analyze this waveguide. Number of Fourier series terms in 
variational method was limited to = 7 and Ny = 23 but can get more 
accurate result by going beyond this limits. 
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Figure 4: Guided and radiation TM-like modes of the crystal of figure 1 
calculated by the variational method and compared with the results of FDTD 
method 

4 Efficiency comparison 

It is clear that the proposed method is much faster than the conventional 
3D PWE method. The reason is that in our method we have to calculate a 
matrix which its dimensions are of the order where N is the number 

of Fourier series terms. While in 3D PWE method the matrix dimensions 
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Figure 5: TM-like magnetic field distribution at the high-symmetry point 
in the mid-plane of the crystal of hgure 1, obtained from (a) variational 
and (b) 3D PWE methods and ^ component of the magnetic held in a vertical 
plane to the crystal obtained from (c) variational and (d) 3D PWE methods. 


are of the order 0{N^). To illustrate this fact quantitatively we have plotted 
in hgure 7 frequencies of the hrst TE-like and fourth TM-like modes at high- 
symmetry points M and K respectively of the introduced crystal in hgure 1 
versus the time it takes to calculate them. These calculations have been per¬ 
formed with a personal computer equipped with a 32-bit Intel® core"'"^2 Duo 
CPU. The frequencies are calculated by variational and 3D PWE methods. 
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Figure 6: Dispersion diagram of the PC slab waveguide shown in the upper 
inset, calculated by both FDTD and variational methods. This waveguide 
is built by Filing one row of air columns in the crystal of hgure 1 in F — K 
direction. The lower inset shows the unit cell used for analysis. Magenta 
regions represent extended modes inside the crystal. 


Results with minimum calculation time are related to = Ny = Nz — 1 = ?> 
limitations on terms of Fourier series and ones with maximum calculation 
time are related to = iVy = — 1 = 8. It can be seen that the varia¬ 

tional method converges to the hnal value at least one order of magnitude 
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faster than PWE method. 



Figure 7: Convergence time comparison between variational and PWE 

methods. Frequencies of the hrst TE-like and fourth TM-like modes at high- 
symmetry points M and K for the crystal of hgure 1 have been plotted versus 
the time it takes to calculate them. A personal computer equipped with a 
32-bit Intel®core™2 Duo CPU has been used to calculate these frequencies. 


5 Conclusions 

We presented a new method for fast analysis of PC slabs. In this method 
which is similar to PWE method instead of creating an artihcial periodicity 
along the normal component of the slab and then approximating the held 
dependency on this component by a Fourier series, we adopted held distri¬ 
bution in a dielectric slab waveguide to estimate held in PC slabs. Results of 
this fast method are in good agreements with that of other standard meth¬ 
ods. We also compared the convergence time of our method to the hnal 
result with that of 3D PWE method. Our method shows at least one order 
of magnitude faster convergence. 
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